A SECOND LOOK AT THE GAUSSIAN SEMICLASSICAL SOLITON 
ENSEMBLE FOR THE FOCUSING NONLINEAR SCHRODINGER EQUATION 



S3 



LONG LEE AND GREGORY D. LYNG 

Abstract. We present the results of a numerical experiment inspired by the semiclassical (zero- 
dispersion) limit of the focusing nonlinear Schrodinger (NLS) equation. In particular, we focus 
£\j ' on the Gaussian semiclassical soliton ensemble, a family of exact multisoliton solutions obtained 

by repeatedly solving the initial-value problem for a particular sequence of initial data. The se- 
quence of data is generated by adding an asymptotically vanishing sequence of perturbations to 
pure Gaussian initial data. These perturbations are obtained by applying the inverse-scattering 
transform to formal WKB approximations of eigenvalues of the associated spectral problem with 
Q ■ a Gaussian potential. Recent results [Lee, Lyng, & Vankova, Physica D 241 (2012) 1767-1781] 

suggest that, remarkably, these perturbations — interlaced as they are with the integrable structure 
of the equation — do not excite the acute modulational instabilities that are known to be present 
in the semiclassical regime. Here, we provide additional evidence to support the claim that these 
WKB- induced perturbations indeed have a very special structure. In particular, as a control experi- 
ment, we examine the evolution from a family of initial data created by an asymptotically vanishing 
^0 ' family of analytic perturbations which are qualitatively indistinguishable from the WKB-induced 

perturbations that generate the Gaussian semiclassical soliton ensemble. We then compare this 
evolution to the (numerically computed) true evolution of the Gaussian and also to the evolution 
of the corresponding members of the semiclassical soliton ensemble. Our results both highlight the 
exceptional nature of the WKB-induced perturbations used to generate the semiclassical soliton 
ensemble and provide new insight into the sensitivity properties of the semiclassical limit problem 
for the focusing NLS equation. 

oo 

00 

1. Introduction 

1.1. Focusing nonlinear Schrodinger equation and modulational instability. We consider 
the initial-value problem for the semiclassically scaled focusing nonlinear Schrodinger (NLS) equa- 
tion 

^40 + 1^0, (1 ,a) 
il>(x,0)=il>o(x). (Lib) 

Here, the dimensionless parameter < e <C 1 measures the ratio of dispersion to nonlinearity, and 
we suppose, to start, that the initial data is given in amplitude-phase form as ipo(x) = Ao(x)e lS °^' e . 
Typical assumptions about the smooth, real-valued functions Aq and Sq are that Aq > and that 
both Aq and S' Q decay rapidly to zero as |x| — > oo. The semiclassical or zero-dispersion limit 
problem is, for fixed initial data ipQ, to solve (jl.ip for each sufficiently small e > and to describe 
the limiting behavior of the resulting family of solutions, denoted ift(x,t;e) and indexed by e, as 
e|0. 

This problem has garnered substantial interest in recent years; for more details and discussion 
we refer the reader to [7J H5j [20j ETJ [26j E3 EH] and the references therein. The mathematical 
interest in this problem is due, in large part, to the phenomenon of modulational instability. To 
describe this phenomenon, we recall that it is a simple matter to verify that a plane wave of the 
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form ip(x,t) = exp(i(kx — ut)) is a solution of (jl.laj) provided that the dispersion relation 



is satisfied. But, a linear stability analysis, perturbing the amplitude and phase slightly, shows 
that, in the limit of vanishing dispersion, all wave numbers are destabilized [24]. As we shall see 
shortly, this feature of equation (jl.la[) has a profound effect on the analysis of the semiclassical 
limit. 

Note 1. The linear stability calculation shows that (jl.lap tends to destabilize periodic wavetrains 
and that this effect becomes more pronounced as e tends to zero. This instability is sometimes 
referred to as the Benjamin-Feir instability [3j [2l] after the discoverers of this phenomenon in the 
setting of water waves. It is also known to occur, for example, in Langmuir waves in plasma [9]. 

Like other examples of zero-dispersion limits in the the theory of integrable nonlinear wave equa- 
tions, many aspects of the analysis semiclassical limit of (jl.lap are best understood in the context 
of the pathbreaking work of Lax & Levermore on the zero-dispersion limit of the Korteweg-de 
Vries (KdV) equation [19]. In Lax & Levermore's analysis, a central role was played by the the 
Whitham (or modulation) equations; these equations describe the local evolution of the large-scale 
structures (the macrostructure) in the highly-oscillatory solutions. Indeed, after an elaborate anal- 
ysis based on the integrability of the KdV equation, Lax & Levermore showed that the weak limits 
of the KdV equation satisfied the Whitham equations (derived independently by Flaschka, Forest, 
& McLaughlin [11]) which describe modulated multiphase KdV wavea^. Importantly, in the case 
of the KdV equation, the Whitham equations are partial differential equations of hyperbolic type. 
Hyperbolicity and its consequent local well-posedness was used by Lax & Levermore in an essential 
way. Namely, they used the local well-posedness to account for a sequence of asymptotically negli- 
gible perturbations they applied to the initial data. These particular perturbations were introduced 
to make the solution of the initial-value problem for the KdV equation by the inverse-scattering 
method especially transparent. 

Given the success of Lax & Levermore's analysis and its integral use of the Whitham equations, 
it is natural to expect that the Whitham systems for (II. lap will be similarly important. In their 
simplest form, one can derive the Whitham equations for (jl.lap by the following intuitive procedure. 
First, we assume that for some positive time (independent of e), the solution can be written in 
amplitude-phase form (like the data) as 



Then, upon making the standard definitions p = A 2 ,p = A 2 d x S, we obtain exactly — that is, with 
no approximation — from (jl.lap the system 



The relevant e-independent initial data is p(x,0) = Aq(x) 2 and p(x,0) = Aq(x) 2 S' (x) . The 
Whitham equations for (jl.lap are those obtained from the system (|1.4p by neglecting the 0(e 2 ) 
term on the right-hand side of (jl.4bp . Clearly, this procedure has some appeal; it eliminates the 
dependence on e in the initial data, and itprovides, at a formal level, what appears to be a reason- 
able model for the semiclassical dynamicaj. However, inspection of the resulting first-order system 

lr The calculations of Flaschka et al. generalized earlier single-phase calculations of Whitham. 

2 That is, the elliptic equations at least have terms that might balance. Obviously, simply setting e to zero makes 
little sense in the original formulation of the problem 

2 






(1.3) 




(1.4b) 



(1.4a) 



of nonlinear partial differential equations reveals immediately that it is of elliptic type and thus is 
typically ill-posed as an initial-value problem. 

Note 2. Given the aforementioned instability of plane-wave solutions in the limit of vanishing e, this 
makes perfectly good sense; the geometric optics ansatz (jl.3p leads to a description of the evolution 
in terms of modulated plane waves, and the aforementioned linear stability analysis suggests that 
plane waves are unstable. 

The above calculations immediately raise questions about the semiclassical (e \. 0) limit for (jl.ip . 
Indeed, it is evident in early attempts to come to terms with this problem that the presence of 
modulational instability led to a great deal of uncertainty about this problem [12J. Early numerical 
simulations [6j Ej showing a disordered region led some to suggest that some kind of statistical 
approach (e.g., averaging as in the study of turbulence) might be required to reveal any kind of 
structure in the apparently chaotic behavior [5J. However, Miller & Kamvissis [23], using a novel 
numerical calculation based on the integrable structure of (jl.lap . presented overwhelming numer- 
ical evidence, at least in the important special case ipo( x ) = 2sechx that, despite modulational 
instability, there is remarkable structure in the limit. Indeed, among the principal results of Kamvis- 
sis, McLaughlin, & Miller [15] is that, for certain analytic data — so that a solution exists by the 
Cauchy-Kovalevskaya theorem, the elliptic Whitham equations indeed give a correct description of 
the macrostructure in the semiclassical limit. 

1.2. Analyticity & Sensitivity in the semiclassical limit. Despite substantial recent progress 
(e.g., |15j 126]). a completely satisfactory theory of the semiclassical limit of (jl.lap remains out of 
reach. To date, the two principal approaches have been based on (i) real, bell-shaped initial data — 
the approach of Kamvissis et al. via semiclassical soliton ensembles described in more detail below; 
and (ii) a particular one-parameter family of data (potentials) for which Tovbis & Venakides [25] 
have been able to treat the associated non-self-adjoint spectral problem (see (11. 5p below) exactly. 
Notably, the latter case features a nonzero phase (So / 0). In either case, given the ellipticity 
of the modulation equations, the importance of analyticity in the data has long been recognized 
in the context of this problem [2]. For example, the numerical experiments of Bronski & Kutz [4] 
demonstrated the extreme sensitivity of the limiting structure to rough perturbations of the data. In 
particular, they found that even a small amount of nonanalyticity (a kink at the origin) completely 
destroyed the regular structure of the peaks and valleys in the post-wave-breaking solution. Building 
on this example, Clarke & Miller [9] showed similar sensitivity even to perturbations of class C 2 (R) 
and to analytic perturbations with singularities off the line in the complex x plane. They found 
that, in the case of the C 2 (M) perturbations, the primary caustic (first wave breaking) was pulled 
back to t = 0, a dramatic demonstration of the sensitivity of the semiclassical limit to the analyticity 
in the data. Given the role of (jl.lap as an important model equation for physical applications, the 
apparent necessity of analyticity in the data is disquieting. Certainly, from the point of view of 
applications, one would like to be able to relax the requirement that the data be at every point 
expandable in power series and to allow for more realistic, rougher data. 

Here we push this line of inquiry in a different direction. In this letter, we consider a family of 
analytic perturbations of real-analytic data designed to mimic the features of perturbations that 
naturally arise in the analysis of the semiclassical limit by soliton ensembles. Notably, Lee, Lyng, 
& Vankova [20] found a remarkable kind of stability in the semiclassical limit in the presence of 
these perturbations. We start by describing the origin of these perturbations; they are intimately 
intertwined with the integrable structure of (jl.lap . and so we begin with a brief description of the 
inverse-scattering transform for (jl.ip . 

1.3. Inverse-Scattering Transform & Semiclassical Soliton Ensembles. For a variety of 
interrelated reasons (including, e.g., the ellipticity of the Whitham system), the semiclassical limit 
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problem for (jl.ip resists any kind of straightforward application of the techniques used by Lax & 
Levermore for the analogous problem for the KdV equation. However, it is also true that the lion's 
share of our current understanding of the small-e behavior of (jl.ip is based on the integrability 
of (jl.lap and the fact that (in principle at least) the initial-value problem may be solved by the 
inverse-scattering transform (1ST). Famously, Zakharov & Shabat [29] showed that (jl.lap has a 
Lax pair, and the 1ST reduces the solution of the nonlinear problem (jl.ip to the solution of the 
(presumably simpler) linear problems of the Lax pair. In the case of (jl.ip . we recall that the 
initial step in this process involves an analysis of the non-self-adjoint Zakharov-Shabat eigenvalue 
problem 

_d 

da 

In (jl.5p . w = {w\,W2) t is an auxiliary function, and A G C is a spectral parameter. Of particular 
interest are those values of A G C for which (|1.5|) has a solution in L 2 (IR). 

In this letter, we restrict our attention to real, bell-shaped initial data. That is, we work in the 
framework of Kamvissis et al. [15] , and we restrict our attention to initial data of the form 

^o(z) = A {x) , (1.6) 

where A$ : M. — > (0, A] is even, bell-shaped, and real analytic. More precisely, Aq is assumed to 
(i) decay rapidly at ±oo; (ii) be an even function, i.e., Aq(x) = Aq(— x) for all (hi) have a 

single genuine maximum at x = 0, i.e., Aq(0) = A, A' Q {Q) = 0, j4/ '(0) < 0; and (iv) be real-analytic. 
For each e > and for = Aq as described above, Klaus & Shaw |17j PT8] have shown that the 
discrete spectrum of (jl.5p is confined to the imaginary axis. They also obtain a lower bound on 
the number of eigenvalues in terms of the size of || A^Il^r)! this bound suggests that the number 
of eigenvalues iV satisfies N^ 1 ~ e. 

Beyond this, it has long been known that in the e ^ limit, the eigenvalue problem (jl.5p can be 
written as the eigenvalue problem for a semiclassical (linear) Schrodinger operator with a (formally) 
small 0(e 2 ) eigenvalue-dependent correction. Indeed, this feature of (|1.5p was noted by Zakharov & 
Shabat [29] in their original paper and later exploited by Ercolani et al. [10] in their analysis. Then, 
neglecting the presumably small correction, one finds that a formal WKB method applied to (|1.5p 
indicates that the reflection coefficient is transcendentally small and that a quantization condition 
of Bohr-Sommerfeld type identifies the locations on the imaginary axis of the, necessarily simple, 
eigenvalues. Since, in general, the true location of the discrete spectrum of (jl.5p is not known, it 
is natural to use the formal WKB values for the discrete spectrum in its place. Similarly, since 
the reflection coefficient is also not typically known but is expected to be small beyond all orders, 
it seems reasonable to replace the true (unknown) reflection coefficient with zero. For each small 
e > this process, neglecting the reflection coefficient and using the WKB eigenvalues induces a 
perturbation in the initial data. That is, the true initial data V'o is replaced with some other initial 
condition ip^ — See Figure [T] below — which depends on e and for which the WKB spectral data 
is the true spectral data. Because reflection is neglected, each solution of (jl.lap with initial data 
ipQ is an iV-soliton with N ~ e _1 . The collection of all these exact multisoliton solutions of (jl.lap 
(with N — > oo and e i 0) is called the semiclassical soliton ensemble (SSE) associated with Aq. We 
write 

and we call q( e > the perturbation induced by the WKB analysis of (jl.5p . For more details, see 

[El 021 [20]. 

This process — the approximation of the initial data by a sequence of reflectionless potentials — 
was exactly that used by Lax &: Levermore |19] in their analysis of the zero-dispersion limit of the 
KdV equation. Following in their footsteps, Kamvissis et al. [15] used the same device as the first 
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step of their analysis. In related work, Miller [22] has shown that ipo and ip are asymptotically 
pointwise close as e 1 0. The important distinction between the two cases is precisely the contrast 
of hyperbolicity versus ellipticity in the corresponding Whitham systems. That is, while Lax & 
Levermore were able to appeal to hyperbolicity to absorb the effect of the vanishing perturbation 
to the data, no such maneuver is available for the semclassical focusing NLS equation. Indeed, 
in the case of the focusing NLS equation, it is necessary to analyze the competition between the 
modulational instabilities whose growth rates become arbitrarily large as e J, and the sizes of the 
perturbations which vanish in the semiclassical limit. 

We omit any further discussion of the inverse-scattering transform. The key point is that, 
once the scattering data (eigenvalues, reflection, proportionality constants^) associated with Aq 
is determined, the other half of the Lax pair enters, and one finds that the scattering data have 
a simple, explicit evolution in time. Thus, finding the solution ip{x, t; e) for t > amounts to 
recovering the new "potential" from the time-evolved scattering data. That is, it is an inverse- 
scattering problem associated with (jl.5p . A very brief outline of this process can be found in the 
Appendix of [20] . 

Note 3 (Notation) . For initial data ipQ , we denote by i/Jq the data induced by the modified scattering 
data as described above. Similarly, we denote the value of the solution with data V>o^ a t the point 
(x,t) by %l)( e \x,t) while we denote by ip(x,t;e) the value of the true solution of (jl.ip with the 
corresponding value of e. Thus, the appearance of a superscript "(e)" always signals the presence 
of an element of the SSE. Later, we shall use a superscript "e" (without parentheses) to denote an 
e-dependent family of perturbations and solutions not associated with the special WKB-induced 
modification of data/solution. Similarly, we shall following the same convention for the density and 
denote by p (or p^ or p e as the case may be) the square modulus of the solution if) (or tj)( € ' or ip e ). 
In most of the graphs in this letter, the vertical axis measures the size of p (or p^ or p e ). 

Recently, Lee et al. [20] computed (by a numerical implementation of the inverse-scattering trans- 
form) many members of the Gaussian SSE, the SSE associated with initial data Aq = exp(— x 2 ). 
These computations required high-precision knowledge of the WKB eigenvalues, a quite delicate 
exercise in root finding, and a numerical implementation of the the 1ST; similar calculations of this 
latter step can be found in |21| I23j. Lee et al. also computed by a finite difference scheme (see §1.41 
below for details about the scheme) the true evolution of the Gaussian initial data Aq = exp(— x 2 ). 
Their numerical experiments showed that the rate of convergence of the perturbations at t = 

|||V>o + g (£) | 2 -^o| 2 || 2 = 0(e) ase|0. (1.7) 

was propagated to small times t > 0: 

\\p^\;t)-p(;t;e)\\ 2 = 0(e) as e | . (1.8) 

(The norm || • || 2 is defined below in (12. 2p .) Thus, despite the presence of modulational instability, 
and the extreme sensitivity displayed by (jl.ip in the presence of rough perturbations to the initial 
data [US], the special perturbations seem not to excite the modulational instability. Panels 
(a) and (b) in Figure [7] show a representative calculation of this type; see also Figure [2] for a visual 
display of the size of the 2-norm difference \\p^ — p\\2 versus e. The solution ip( e ) evolving from the 
perturbed data tracks the true solution •; e) remarkably well, even past the primary caustic. 

In this letter, we construct a vanishing family of real analytic perturbations satisfying (jl.7p . 
Our construction has a dual motivation. First, the constructed perturbations are qualitatively 
indistinguishable from the family ■ We thus view the computations in this letter as a control 



The proportionality constant associated with each eigenvalue carries the requisite information about the associated 
eigenfunction. 
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experiment to complement the computations in [20]. Indeed, we find that the preservation of the 
rate of convergence of the perturbations to zero appears to be far from a generic property of the 
evolution under (jl.lap . This bolsters the claims of |20| . and it highlights the special structure of 
the WKB-induced perturbations. In particular, it gives additional evidence to support the claim 
that the WKB approximation of the spectrum of (|1.5p is valid. Our secondary motivation comes 
from the fact that there are now the beginnings of a robust theory for the semiclassical limit of 
the focusing NLS equation with analytic data; we hope this letter pushes the line of inquiry in the 
direction originally posed as a somewhat open-ended question by Kamvissis et al. [15]. That is, 
Kamvissis et al. wondered about the sensitivity of the limit within the class of analytic data (or at 
least within some class of admissible analytic perturbations). For example, they wondered about 
the stability of the limit under small (in L 2 (M), say) analytic perturbations. We believe that our 
experiment in this letter represents a first, preliminary step in this direction. The result of [20] 
suggests that there is at least one family of analytic perturbations under which the limiting behavior 
is quite robust. The determination of a class of "admissible perturbations" and a description of 
this class (even for analytic perturbations) for the semiclassical limit problem remains a substantial 
open problem. 

1.4. Numerical Methods. The construction of accurate numerical methods for approximating 
the solution of (11. ip when e is small is a notoriously difficult problem. We refer the interested 
reader to the recent survey of Jin et al. [13J for a comprehensive discussion of the various challenges 
that occur even in the linear case; see also [U [U EJ [20] for additional details about numerical 
methods for (jl.ip when e is small. Among the challenges associated with computing high-quality 
approximate solutions of (jl.ip . we recall that the solution of (jl.ip is expected to develop small 
structures of size O(e) in time and space (the microstructure). Thus, when e is small, one would 
expect that a very small spatial mesh size Ax and time step At need to be used to completely 
resolve the microstructure. However, the high-wavenumber modes (short wavelengths) introduced 
by the small mesh size could be exacerbated by modulational instability and quickly destroy the 
fidelity of the approximation. Indeed, as described in [20], the interplay between these aspects of 
the computation can be subtle. 

For all of the calculations in this letter, we use a second-order implicit finite-difference algorithm 
developed in [20] and shown there to be a suitable numerical method for solving the semiclassical 
focusing nonlinear Schrodinger equation (jl.ip in the small-e regime. The algorithm adopts the 
completely integrable spatial discretization proposed by Ablowitz & Ladik [TJ. Then, the midpoint 
time integrator is applied to advance the system in time; a simple fixed-point iteration (FPI) is 
used at each time step to obtain the update for the next time step. The initial guess for the FPI 
procedure is the solution obtained by solving the ODE by using a Crank-Nicolson type of method 
(trapezoidal method for the linear component and the Euler method for the nonlinear one). More 
details about the algorithm, including its validation and assessment, can be found in the original 
paper [20] . 

2. Gaussian SSE 

2.1. The Structure of the Gaussian SSE. As in [20], we now restrict our attention to the case 
that 

V>o(x) = A (x) = e~ x2 . (2.1) 

As described above, Lee et al. |20] recently investigated the effect of the perturbations that give rise 
to the SSE associated with ipQ. First, we note that Figure Q] shows representative reconstructions of 
the square modulus of rp^ for two different values of e. As noted above, these reconstructions were 
done using a numerical implementation of the inverse-scattering transform as in |23l 1201 \2l\ . As in 

those previous cases, the multisoliton solutions associated with data ip^ are completely specified 
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by the eigenvalues and the proportionality constants; in such cases, because reflection is neglected, 
the reconstruction (1ST) reduces to the solution of a linear system. The linear system must be 
solved independently for each x and t; its size and condition number grow as e \, 0. 

Note 4. Evidently, despite the fact that the WKB eigenvalues used in the construction of ip^ are 
purely imaginary, the modified data shown in Figure [Q is not of Klaus-Shaw type [TT| [18] . 



Modified initial data with SSE, £=1/(5V(n)) 



Modified initial data with SSE, e=1/(15\'(7r)) 





-2.5 -2 -1.5 -1 -0.5 0.5 1 1.5 2 2.5 



-2.5 -2 -1.5 -1 -0.5 0.5 1 1.5 2 2.5 



(a) " (b) 

Figure 1. Modified initial data (t = 0) for SSE: (a) e = l/(N^),N = 5; (b) 
e = l/(NT,/Tr),N = 15. The solid lines are the modified data, and the dashed lines 
are the unmodified Gaussian data i/jq. 

Figure [2] shows the result of the comprehensive comparison between ip^ and t/j carried out in 
[20] . That is, the figure shows that the error for t = 0.0,0.1, . . . ,0.5 between the members of the 
SSE and numerical approximations of the true solution with the corresponding values of e. In the 
figure, the error is computed using the norm on computed quantities defined by 



\U\\ 



1 M 



(2.2) 



where M is the number of spatial grid points considered and Ui is the U value at the i th spatial 
grid point. Thus, the vertical axis in Figure [2] is 



E = \\\^\-,t)[ A - m;t;e N ){% = \\p^>(.,t)-p(;t;e N )\\ 2 , 
where f// ejv ) is a member of the SSE (an N soliton), and ip represents the evolution of (jl.ip with 

2 

■00 = e _x and e = e^. The structure of the error, particularly, the consistent decay at an 0(e) 
rate, is striking about this experiment. 

2.2. Finite-difference evolution for perturbed initial data ip^ . Before we perform the prin- 
cipal experiment of this letter, we take the oscillatory functions as initial data, and we verify 
that the finite-difference scheme described in §1.41 produces good approximations to the solution 
of the focusing NLS equation with such data for the ranges of e and t that we consider in this 
letter. To this end, we will compare the finite-difference approximation evolved from initial data 
i/jq^ with the 1ST calculation (described above) at a given final time, and we will show that they 
are indeed in good agreement. To reiterate, here we are computing the same object in two different 
and completely independent ways. As noted above, the 1ST calculation is performed independently 
for each x and t, whence numerical errors do not propagate. Figure shows that both methods 
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Figure 2. The data points show the 2-norm error E between \ift( eN '(-, i)| 2 and 
|^(-, £at)| 2 with 6 = l/(y/TrN) against iV = 1,2,..., 20 for time slices t = 
0,0.1,0.2,0.3,0.4,0.5. The curves are of the form E = C ■ N a where the constants 
C and a are determined by least squares. The values of a in the legend show an 
0(N~ 1 ) = O(e) rate of convergence at t = and for later times too (even past the 
primary caustic). This figure is taken from |20j . 

(finite differences and 1ST) of obtaining iPq\-, 0.5) give answers which are in strong agreement. The 
small parameter takes the values e = 1/(5-^/77) ~ 0.1128 and e = 1/(15-^/77) ~ 0.01175 in (a) and 
(b), respectively. Indeed, the 2-norm difference of p between the two calculations is 1.5766E-05 for 
Figure Ela), and is 2.2801E-3 for Figure [3jb). The differences are measured for —1 < x < 1. The 
good agreement shown in Figure [3] strongly supports the claim that the evolution of the perturbed 
Gaussian initial data produced by our finite difference scheme is not some numerical artifact. 



indexed by e > 0. We shall use p e as an e-dependent perturbation of the initial data f/'oj we 
refer to any member of this family of functions as a "cosine perturbation." Figure [4] shows the 
modification to the initial data using the cosine perturbation p e in equation (13. ip for e = 1/(5^71") 
and e = 1/(15-^/7?), respectively. We note that the values 0.3 and 0.54 that appear in (13. ip are 
chosen so that the amplitude and frequency of the cosine perturbations p e track their counterparts 
g( e ); compare Figure H] and Figure [TJ From a qualitative point of view, these data look very similar 
to the modified data for the SSE in Figure [TJ Indeed, without the Gaussian template, they are 
virtually indistinguishable by eye. To more clearly illustrate the nearness of the perturbations, 
Figure [5] shows a plot of the difference 



3. Experiment 



Consider the family of functions 
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SSE modified initial data, 6=1/(5 Vjc), time = 0.5 



SSE modified initial data, e=1/(15 Vsc), time = 0.5 



SSE (1ST) _ 

FD evolution 



SSE (1ST) 

FD evolution 
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-1 -0.8 -0.6 -0.4 -0.2 0.2 0.4 0.6 0.8 



(a) " (b) 

Figure 3. Finite-difference evolution from the modified initial data for SSE (Figure 
QJ. (a) e = l/(Ny/Tr),N = 5, at t = 0.5, (b) e = 1/(Nt/tF), N = 15, at t = 0.5. The 
dashed lines are the finite-difference evolution from the modified initial data, and 
the solid lines are the corresponding Gaussian SSE at t = 0.5. The difference of p 
for (a) in the 2-norm is 1.5766E-05 measured for —1 < x < 1. Similarly the 2-norm 
difference for (b) is 2.2801E-3. 



for x 6 [—2.5,2.5] again for e = l/(5\/7r) and e = l/(15^/^v). A straightforward calculation shows 
that if < e < 1, then 

/oo 
|p e (x)| 2 dx < e 2 x (0.07) (1 + e" 1 ' 8 ) , (3.2) 
-oo 

and thus ||p e ||.L2( R ) = O(e) as e J, 0. Observe also that the cosine perturbation p e is even (like 
q( e t) and oscillatory with frequency depending on e (like q^)- Moreover, for each fixed e, p £ is the 
composition of real analytic functions hence is real analytic. 

Now, the basic question is whether or not the evolution from the family of initial data 

V>o = ipo + P € 

has structure as e tends to zero that matches the remarkable structure, shown in Figure [21 of 
the solutions obtained from the data Vrj^ = V'o + ■ The answer is a definitive no. Figure [6] 



cosine perturbed initial data,f=1/(5v{ir)) 



cosine perturbed initial data,E=1/(15^(jr)) 




-2.5 -2 -1.5 -1 -0.5 0.5 1 1.5 2 2.5 



(a) 




2.5 -2 -1.5 -1 -0.5 0.5 1 1.5 2 2.5 



(b) 



Figure 4. Modified initial data using the cosine perturbation (Eq. (|3. lfl ) : (a) 
e = (b) e = l/ilhyfw). The solid lines are the modified data, and the 

dashed lines are the unmodified Gaussian. 
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Figure 5. The difference in the magnitude of p between the SSE and the cosine 
modified initial data (Figures [L] and H|) : (a) e = l/(5y^); (b) e = l/(15i/7r). 



shows plots of the two densities p e and p at t = 0.3 corresponding to the initial data including 
the perturbation (I3.1j) and to the pure Gaussian data (12. lj) . In the figure, the small parameter e 
takes the representative values where N = 5,10,15 and 20, respectively. The chosen time 

slice, t = 0.3, is before the upper bound on the breaking time found in [2D] for the Gaussian 
initial data. Evidently, Figure [6] shows that, in marked contrast of the case of p( € \ the solution p e 
evolving from vfjQ + p € develops small-scale oscillations well in advance of the breaking time for the 
unperturbed problem. Indeed, the behavior displayed in Figure [6] suggests that, unlike the SSE 
case in |20j . systematically computing the 2-norm difference between p and p e for a sequence of 
decreasing values of e might not reveal any meaningful information. 

To further illustrate the behavior of the cosine perturbation, we evolve, using the finite-difference 
scheme, the SSE initial data (Figure QJb)) and the cosine perturbed data ipQ (Figured] (b)) to 
t = 1, and we compare the results with the evolution of the Gaussian initial data ipQ. Figure [7] shows 
the result for e = l/(15-y/7r). Panel (a) shows the square modulus of the initial data ipo( x ) an d the 
corresponding surface plot of \ip(x, t; e)| 2 over the unit square [0, 1] x [0, 1] in the xi-plane. Similarly, 
panels (b) and (c) show the data and solution corresponding to the WKB-induced perturbation 
q( e > and the artificial, but strikingly similar, perturbation p e , respectively. The figure clearly shows 
that the evolutions of the smooth Gaussian initial data and the SSE initial data are virtually 
indistinguishable, while the evolution of the cosine perturbed initial data is drastically different 
from the other two cases. Indeed, the early onset of oscillations in advance of the breaking time for 
the smooth Gaussian is clearly visible; this behavior is somewhat reminiscent of the phenomenon of 
"beards" described by Clarke & Miller [9j in an experiment which featured initial data with points 
of non-analyticity near the real x axis. 

4. Discussion 

The principal conclusion of our experiment is that we see no structure of decay in the error 

\\p*(-,t)-p(-,t-e)\\ 2 ; 

that is, we see no pattern that appears to be in any way close to that shown in Figure[2] We interpret 
this as evidence, bolstering the claim of |20] . that the observed O(e) decay in the case of the SSE is 

truly due to the special structure embedded in ip^ and not due to its coarser features (i.e., those 
present in i/jq as well). To say this another way, one suspects that the WKB approximations to 

the eigenvalues of (II. 5h used to construct ip^' are indeed very good approximations to the true 
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Figure 6. Comparison of the amplitude of p at t = 0.3 between the cosine 
perturbed initial data and the Gaussian initial data, (a) e = l/(5- v /7r); (b) 
e = 1/(1(V5F); (c) e = 1/(150F); (d) e = 1/(20 Vi). 



(unknown) eigenvalues associated with the Gaussian. Indeed, we see a rigorous proof of this claim 
as a very important open problem. Of course, similar reasoning leads one to believe, based on the 
computed behavior of the solution ip e , that the scattering data associated with the artificial data 
tpQ is not so close to that of the pure Gaussian. The issue, of course, is that at present it is not 
at all clear in general how to extract all of the necessary information about the scattering data for 
a given potential. In the current setting, one would also like to understand the sensitivity of the 
spectrum to perturbations of the potential. 

Conventional wisdom, based on formal analysis of (|1.5|) . seems to hold that, compared to the 
general case of a complex-valued potential, the structure of the spectrum of fjl .5[) is more stable in 
case of a real potential tpQ. For example, the aforementioned reduction of the eigenvalue problem 
to that of a perturbed Schrodinger operator and the behavior of the famous unpublished "shadow 
bound" of Deift, Venakides, & Zhou together give credence to this intuition. Moreover, numerical 
calculations of the spectrum by Bronski [5] showed a remarkable stability in the real case even in 
response to a non-analytic perturbation. However, we note that Bronski considered the potential 
Aq(x) = sechx together with perturbation r e (x) = eexp(— \x\), and the combination 



4> e = A + r e 
n 



True Data 



Gaussian, E = 1/(15Vit) 



smooth Gaussian 



(a) 



(b) 



(c) 



-2.5 -2 -1 .5 -1 -0.5 0.5 1 1 .5 2 2.5 



Modified initial data with SSE, e=1/(15\'(ti)) 




-SSE initial 

- smooth Gaussian 



-2.5 -2 -1.5 -1 -0.5 0.5 1 1.5 2 2.5 



perturbed initial daia.e=1/(15V(jt)) 




-1.5 -1 -0.5 




t (time] 



SSE, e = 1/(15V:t) 




t (time] 



e perturbation, e = 1 /(1 5Vjt) 




t (time] 



Figure 7. LEFT: Initial data (a) Pure Gaussian, i(j = exp(— x 2 ); (b) SSE at t = 0, 

Wi^ = "00 + (c) Cosine perturbation, V'o = RIGHT: The time-space 

plot for (a) Gaussian initial data, \ip(x, t; e)| 2 ; (b) the SSE modified initial data, 
\ip( e \x , t)\ 2 ; and (c) the square modulus of the solution evolved from the cosine 
perturbed initial data, \ip e (x, t)\ 2 . In this figure e = 1/(15^/7?) « 0.01175 for all 
three cases. 



belongs to the class of Klaus-Shaw potentials |17} [18jj. One might reasonably wonder, in the 
presence of the experiment presented here, whether or not the stability phenomenon observed by 
Bronski is more a feature of the fact that the perturbed potential is of Klaus-Shaw type than its 



4 The key requirements utilized by Klaus & Shaw to obtain confinement of the spectrum to the imaginary axis are 
that t/jq be a non-negative, bounded, piecewise-smooth, square-integrable function that is nondecreasing for x < 
and nonincreasing for x > 0. 
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reality. We emphasize that for small e the initial data corresponding to the cosine perturbation, 
ipQ is clearly not in the Klaus-Shaw [T7] class of initial data; the critical "concentration" property 
does not hold. Indeed, Klaus & Shaw [16] have shown that their confinement result can fail for 
symmetric "double-lobe" potentials. Perhaps that is the case here. However, given the oscillations 
present in ipQ, the failure of the oft-made assumption of two turning points makes even the formal 
analysis, as in [ID], unclear. Indeed, motivated by the desire to avoid the complexity of such a 
WKB-based analysis, [27] have recently proposed an alternative approach, based on an integral 
transform that can be viewed as a complexified version of the Abel transform, to extracting from 
ipo the leading-order term of its scattering data. 

Finally, we hope that our experiment will serve to stimulate a broader discussion of the sensitivity 
of the semiclassical limit for even in the analytic class. For example, an interesting and closely 
related follow-up problem would be to investigate the continuity properties of one or more of the 
important features of the semiclassical limit (e.g., the t coordinate of the first breaking at x = 0) 
with respect to perturbations in the space Fjc — the space of band-limited function^ — proposed by 
Clarke & Miller |9j. 
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